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Q , In this paper, we propose a novel large deformation diffeomorphic registration algorithm to align 

high angular resolution diffusion images (HARDI) characterized by orientation distribution functions 
^ ■ (ODFs). Our proposed algorithm seeks an optimal diffeomorphism of large deformation between two 

m 

\^ ' ODF fields in a spatial volume domain and at the same time, locally reorients an ODF in a manner such 



that it remains consistent with the surrounding anatomical structure. To this end, we first review the Rie- 



r~^ . mannian manifold of ODFs. We then define the reorientation of an ODF when an affine transformation 

^■"■^ ■ is appUed and subsequently, define the diffeomorphic group action to be applied on the ODF based on 



this reorientation. We incorporate the Riemannian metric of ODFs for quantifying the similarity of two 
HARDI images into a variational problem defined under the large deformation diffeomorphic metric 
mapping (LDDMM) framework. We finally derive the gradient of the cost function in both Riemannian 
spaces of diffeomorphisms and the ODFs, and present its numerical implementation. Both synthetic and 
real brain HARDI data are used to illustrate the performance of our registration algorithm. 
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Orientation distribution function (ODF), diffeomorphic group action on ODF, large deformation 
diffeomorphic metric mapping, ODF reorientation. 

I. Introduction 

Diffusion weighted magnetic resonance imaging (DW-MRI) is a unique in vivo imaging 
technique that allows us to visualize the three-dimensional architecture of neural fiber pathways 
in the human brain. Several techniques may be used to reconstruct the local orientation of 
brain tissue from DW-MRI data. A classical method is known as Diffusion Tensor Imaging 
(DTI) [jll, which characterizes the diffusivity profile of water molecules in brain tissue by a 
single oriented 3D Gaussian probability distribution function (PDF). In DTI, the diffusivity 
profile is often represented mathematically by a symmetric positive definite (SPD) tensor field 
D : M^ — > SPD(3) C M^^^ that measures the extent of diffusion in any direction v G M^ 
as v^Dv. The geometry of SPD(3) is well-studied and several metrics for comparing tensors 
have been proposed [|2l, flSl, dH, flSl. Based on these metrics, statistical tests such as voxel-based 
analysis of diffusion tensors have been developed JH, 0, flU, flU. Before such population studies 
can been carried out, there is a essential need to perform DTI registration, that is, to align tensor 
data across subjects to a standard coordinate space. 

Compared to the classical image registration problem, the registration of DTI fields is more 
complicated since DTI data contains structural information affected by the transformation. Two 
key transformations need to be defined: a transformation to spatially align anatomical structures 
between two brains in a 3D volume domain, and a transformation to align the local diffusivity 
profiles defined at each voxel of two brains. More precisely, a transformation of the image 
domain induces a reorientation of the DTI as the direction of diffusion depends on the coordinate 
system. Thus, for two diffusion tensors Di(a;) and D2(x) at voxel x, it is no longer true that 
Di(x) ~ D2(0(x)) and each tensor must be reoriented in such a way that it remains consis- 
tent with the surrounding anatomical structure. There exist several approaches for reorientation 
that are used in DTI [fTOl . For instance, the Finite Strain (FS) scheme decomposes an affine 
transformation matrix A into A = RS, where R is the rigid rotation and S is the deformation, 
and reorients the tensor D as RJDR^ . An alternative strategy is the Preservation of Principal 
Direction (PPD), in which the reoriented tensor keeps its eigenvalues, yet its principal eigenvector 
vi is transformed as Avi/pvi||. The reader is referred to ifTTIl. [JHll. llBll. [IT4]|. lIBll. [IT6]| and 
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references therein for the existing DTI registration methods. 

While it has been demonstrated that DTI is valuable for studying brain white matter de- 
velopment in children and detecting abnormalities in patients with neuropsychiatric disorders 
and neurodegenerative diseases, a major shortcoming of DTI is that it can only reveal one 
dominant fiber orientation at each location, when between one and two thirds of the voxels 
in the human brain white matter are thought to contain multiple fiber bundles crossing each 
other ifTTl . High angular resolution diffusion imaging (HARDI) ifTSl addresses this well-known 
limitation of DTI. HARDI measures diffusion along n uniformly distributed directions on the 
sphere and can characterize more complex fiber geometries. Several reconstruction techniques 
can be used to characterize diffusion based on the HARDI signals. One class is based on higher- 
order tensors [[T9ll , ll20ll and leverage prior work on DTI. Another method is Q-ball Imaging, 
which uses the Funk-Radon transform to reconstruct an orientation distribution function (ODF). 
The model-free ODF is the angular profile of the diffusion PDF of water molecules and has 
been approximated using different sets of basis functions such as spherical harmonics (SH) 
EH, Ell, E3l, EH, El- Such methods are relatively fast to implement because the ODF 
is computed analytically. By quantitatively comparing fiber orientations retrieved from ODFs 
against histological measurements, Leergaard et al. [26] shows that accurate fiber estimates can 
be obtained from HARDI data, further validating its usage in brain studies. 

Similar to the case of DTI, an open challenge in the analysis of mathematically complex 
HARDI data is registration. Several HARDI registration algorithms have been recently proposed 
under a specific model of local diffusivity. Chiang et al. Il27l proposes an information-theoretic 
approach for fluid registration of ODFs. An inverse-consistent fluid registration algorithm that 
minimizes the symmetrized KuUback-Leibler divergence (sKL) or J-divergence of the two DT 
images [[T6l is first performed and the ODF fields are registered by applying the corresponding 
DTI mapping. The ODFs are reoriented using the PPD method where the principal direction of 
the ODF is determined by principal component analysis. Cheng et al. Il28l takes the approach 
of representing HARDI by Gaussian mixture fields (GMF) and assumes a thin-plate spline 
deformation. The L^ metric of GMFs is minimized, and reorientation is performed on the 
individual Gaussian components, each representing a major fiber direction. Barmpoutis et al. ||29l 
uses a 4th order tensor model and assumes a region-based nonrigid deformation. The rotationally 
invariant Hellinger distance is considered and an affine tensor reorientation, which accounts for 
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rotation, scaling and shearing effects, is applied. Geng et al. Il30l performs a diffeomorphic 
registration is performed with the L^ metric on ODFs represented by spherical harmonics. 
Reorientation is done by altering the SH coefficients in a manner similar to the FS method 
in DTI where only the rotation is extracted and applied. Bloy et al. ll3T| performs alignment of 
ODF fields by using a multi-channel diffeomorphic demons registration algorithm on rotationally 
invariant feature maps and uses the FS scheme in reorientation. Yap et al. [[32ll uses the SH- 
based ODF representation and proposes a hierarchical registration scheme, where descriptors 
are extracted at each level and the alignment is updated by using features extracted from the 
increasing order of the SH representation. Reorientation is done by tilting the gradient directions 
via multiplying with the local affine transform and normalizing. 

Paper Contributions. Unlike a majority of the above-mentioned HARDI registration approaches 
that seek small deformation between two brains, we present a novel registration algorithm for 
HARDI data represented by ODFs under the framework of large deformation diffeomorphic 
metric mapping (LDDMM) such that the deformation of two brains is diffeomorphic (one-to- 
one, smooth, and invertible) and can be in a large scale. Previous studies ll33l suggest that the 
transformation from one brain to another can be really large and therefore small deformation 
models may not be enough. Our proposed algorithm seeks an optimal diffeomorphism of large 
deformation between two ODF fields across a spatial volume domain and at the same time, locally 
reorients an ODF in a manner that it remains consistent with the surrounding anatomical structure. 
We define the reorientation of an ODF when an affine transformation is applied and subsequently, 
define the diffeomorphic group action to be applied on the ODF based on this reorientation. The 
ODF reorientation used in this paper ensures that the transformed ODF remains consistent with 
the surrounding anatomical structure and at the same time, not solely dependent on the rotation. 
Rather, the reorientation takes into account the effects of the affine transformation and ensures 
the volume fraction of fibers oriented toward a small patch must remain the same after the patch 
is transformed. The Riemannian metric for the similarity of ODFs is then incorporated into a 
variational problem in LDDMM. Finally, we derive the gradient of the cost function in both 
Riemannian spaces of diffeomorphisms and the ODFs and present its numerical implementation. 
Even though this paper is based on our previous work ll34l . one major fundamental difference is 
that the gradient derivation in this paper account for orientation differences in the ODFs while 
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||34l does not. We will elaborate how the proposed algorithm outperforms that in ||34ll while 
we discuss the gradient derivation in in §II-E[ Our experiments are shown on synthetic and real 
HARDI brain data in gml 

II. Methods 

A. Review: the Riemannian Manifold of ODFs 

As mentioned in Jl HARDI measurements can be used to reconstruct the ODF, the angular 
profile of the diffusion probability density function (PDF) of water molecules. The ODF is 
actually a PDF defined on a unit sphere S^ and its space is defined as 

P = {p : §2 ^ M+|Vs e §^p(s) > 0; / p{s)ds = 1} . 

The space of p forms a Riemannian manifold, also known as the statistical manifold, which is 
well-known from the field of information geometry [|35l . Rao |[36l introduced the notion of the 
statistical manifold whose elements are probability density functions and composed the Rieman- 
nian structure with the Fisher-Rao metric. Il37l showed that the Fisher-Rao metric is the unique 
intrinsic metric on the statistical manifold P and therefore invariant to re-parameterizations of 
the functions. There are many different parameterizations of PDFs that are equivalent but with 
different forms of the Fisher-Rao metric, leading to the Riemannian operations having different 
computational complexity. In this paper, we choose the square-root representation, which is 
used recently in ODF processing ll38l . Il39l . The square-root representation is one of the most 
efficient representations found to date as the various Riemannian operations such as geodesies, 
exponential maps, logarithm maps are available in closed form. 



The square-root ODF (vODF) is defined as il^is) = a/p(s), where il^{s) is assumed to be 
non-negative to ensure uniqueness. The space of such functions is defined as 

* = {V : S^ ^ M+|Vs e §2^ V(s) > 0; / il)\s)ds = 1}. (1) 

Jses2 

We see that from Eq. ©, the functions t/> lie on the positive orthant of a unit Hilbert sphere, 
a well-studied Riemannian manifold. It can be shown ||40l that the Fisher-Rao metric is simply 
the L^ metric, given as 

(^,-,6)v'. = / li(s)6(s)t/s, (2) 

JsGS2 
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where ^j, ^k £ T^^^ ^e tangent vectors at i/^j. The geodesic distance between any two functions 
ipi, ipj G ^ on a unit Hilbert sphere is the angle 

dist(V'i, i/'j) = II logw,,(V'j)llv, = cos'^itpi, ipj) = cos"^ ( / xpi{s)il)j{s)ds ] , (3) 

where (■, ■) is the normal dot product between points in the sphere under the L^ metric. For the 
sphere, the exponential map has the closed-form formula 

exp^,(0 = cosdl^ll^JV'i + sin(||^||v,J j^, (4) 



where $, G T^,* is a tangent vector at Vi and ||^||^^ = a/(|7I)v^- By restricting |||||^^ G [0, f], 
we ensure that the exponential map is bijective. The logarithm map from ■i/jj to Vi has the 
closed-form formula 



^ = iog^^(^,) = \ ^^;'^^^^^; cos-i(V'. ^,). (5) 

B. Affine Transformation on Square-Root ODFs 

In this section, we discuss the reorientation of the VODF, 'tp{s), when an affine transformation 



A is applied. We denote the transformed VODF as "0(8^) = Aip{s), reflecting the fact that an 
affine transformation induces changes in both the magnitude of ip and the gradient directions of 
s. We will now illustrate how the reorientation is done. 

First of all, we discuss the change in the gradient directions of s. We assume that the change 
of the gradient directions due to affine transformation A is 

- -^"'^ (6) 



\\A-^s\\ ' 

where the transformed gradient directions s^ are normalized back into the unit sphere §^. Notice 
that for s G §^, Eq. Q defines an invertible function of s and therefore, we can find the ODF 
Atl}{s) using the change-of- variable technique of PDF. This will give us the following theorem. 

Theorem II.l. Reorientation of -«/' based on affine transformation of A. Let Ail^{s) be the 

result of an affine transformation A acting on a yODF ip{s). The following analytical equation 
holds true 






^'*w=,/^^^:^ ^{,f^\, (7) 
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where \\-\\ is the norm of a vector 



The ODF reorientation used in this paper ensures that the transformed ODF remains consistent 
with the surrounding anatomical structure and at the same time, not solely dependent on the 
rotation. Rather, by constructing the change-of-variable technique, the reorientation takes into 
account the effects of the affine transformation and ensures the volume fraction of fibers oriented 
toward a small patch must remain the same after the patch is transformed. Figure [T] illustrates 
how A%j){s) varies when A is a rotation, shearing, or scaling and ■j/'(s) contains a single fiber or 



crossing fibers. By construction, Axl){s) fulfills the definition of the yODF. Hence, the similarity 
of A%l>{s) to the square-root ODFs can be quantified in the Riemannian structure given in §II-AI 
for the HARDI registration. 

C. Diffeomorphic Group Action on Square-Root ODF 

We have shown in §II-BI how to reorient ?/? located at a fixed spatial position x in the image 
volume 17 C M'^ through an affine transformation. In this section, we define an action of 
diffeomorphisms (p : CI — )■ Ct on tp, which takes into consideration the reorientation of -0 



as well as the transformation of the spatial volume in 17 C M^. Denote ip{s,x) as the VODF 
with the orientation direction s G §^ located at a; G f7. We define the action of diffeomorphisms 
on ip{s, x) in the form of 

(f) ■ V(S, X) = A^-1(^)'0(S, (j)~^{x)), 

where the local affine transformation A^ at spatial coordinates x is defined as the Jacobian matrix 
of evaluated at x, i.e., A^ = Dx4>- According to Eq. (|7]), the action of diffeomorphisms on 
■0(s, x) can be computed as 



■ V(s,x) 



det [Dip-K 



, -1 



[D^-i(h) s 




(8) 



For the sake of simplicity, we denote (f) ■ ^jj{s, x) as 

(j)-ip{s,x) = Aipo(f)-^{x) , (9) 

where it will be used in the rest of the paper. 
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(a) Single Fiber 

^w ^^ ^^ 

Py 0.1 0.2 0.5 

Shearing ^ ^ ^ ^ ^0 


Sy 0.8 1.2 1.5 

(b) Crossing Fibers 


Rotation 


0^ 15° 30° 45" 


Shearing 
Scaling 


Py 0.1 0-2 0.5 
^V 0.8 1.2 1.5 



Fig. 1: Examples of local affine transformations on ODFs witii a single orientation fiber (panel (a)) and crossing fibers (panel (b)). 
From top to bottom of each panel, three types of affine transformations. A, on the ODFs are demonstrated: a rotation with angle 
9z, where A = [cos^^ — sin^z 0; sin^z cos ^2 0; 001]; a vertical shearing with factor py, where A = [1 0; —py 1 0; 001]; 
and a vertical scaling with factor <;y where yl = [1 0; <iy 0; 1] . 
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Since ■ if){s,x) is in the space of VODF, the Riemannian distance given in §II-AI can be 



directly used to quantify the similarity of (p ■ ^|}{s, x) to other yODFs, which we employ in the 
HARDI registration described in the following section. 

D. Large Deformation Diffeomorphic Metric Mapping for ODFs 

The previous sections equip us with an appropriate representation of the ODF and its diffeo- 
morphic action. Now, we state a variational problem for mapping ODFs from one volume to 
another. We define this problem in the "large deformation" setting of Grenander's group action 
approach for modeling shapes, that is, ODF volumes are modeled by assuming that they can be 
generated from one to another via flows of diffeomorphisms (pt, which are solutions of ordinary 
differential equations (j)t = Vt{(t)t),t E [0, 1], starting from the identity map 00 = Id. They are 
therefore characterized by time-dependent velocity vector fields Vt,t E [0,1]. We define a metric 
distance between a target volume i/'targ and a template volume t/^tcmp as the minimal length of 
curves 0* ■ V'tcmp, t E [0, 1], in a shape space such that, at time t = 1, 0i ■ -01001? = V'targ- Lengths 
of such curves are computed as the integrated norm \\vt\\v of the vector field generating the 
transformation, where Vt E V, where V is a reproducing kernel Hilbert space with kernel ky 
and norm || ■ ||y. 

To ensure solutions are diffeomorphisms, V must be a space of smooth vector fields (4V\. 
Using the duality isometry in Hilbert spaces, one can equivalently express the lengths in terms 
of rrit, interpreted as momentum such that for each u eV, 

{mt, u o 0^)2 = {ky^vt, u)2, (10) 

where we let (m, u)2 denote the L^ inner product between m and u, but also, with a slight abuse, 
the result of the natural pairing between m and v in cases where m is singular (e.g., a measure). 
This identity is classically written as (plrrit = ky^vt, where (p^ is referred to as the puUback 
operation on a vector measure, m^. Using the identity \\vt\\v = (^\/^^t)^t)2 = ("^i5^v"^t)2 and 
the standard fact that energy-minimizing curves coincide with constant-speed length-minimizing 



curves, one can obtain the metric distance between the template and target vODF volumes, 

p(V'tcmp, V'targ), by minimizing /^ {nit, kvmt)2dt such that 0i ■ V'tcmp = V'targ at time t = I. 
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We associate this with the variational problem in the form of 

J{mt) = inf Pl-^temp, V'targ)^ + A / ^3.(01 ■ V'tcmp(s,x),?/>targ(s,x))rfx (11) 

mt:<j>t=kvmt{<f>t), J x^ii 

with Ex as the metric distance between the deformed VODF template, 0i ■ V'tcmp(s, x), and the 

target, ipta.vg{^,x). We use the Riemannian metric given in §II-AI and rewrite Eq. (fTTI) as 

1 p 

2 



J{mt)= inf {mt,kvmt)2dt + X W^og^^p o^-\4^t^rs{x))\\ , dx, 

mt:<f>t=kvmt(^t),Jo J xGil 

0O=Id 

(12) 
where A = Dcpi, the Jacobian of 0i. For the sake of simplicity, we denote 'j/'targ(s, x) as ■j/'targ(a^)- 
Note that since we are dealing with vector fields in M'^, the kernel of V^ is a matrix kernel operator 
in order to get a proper definition. We define this kernel as kyld^xs, where Idsxs is an identity 
matrix, such that ky can be a scalar kernel. In the rest of the paper, we shall refer to this 
LDDMM mapping problem as LDDMM-ODF. 

E. Gradient of J with respect to rrit 

The gradient of J with respect to mj can be computed via studying a variation ml = mt + efht 
on J such that the derivative of J with respect to e is expressed in function of mj. According 
to the general LDDMM framework derived in [|42|. [|43l . we directly give the expression of the 
gradient of J with respect to rrit as 

VJ{mt) = 2mt + Ar/i , (13) 

where 

r]t = V^^E+ [d^A^Vrris)] {r]s + ms)ds , (14) 

where 9^^(/i;yms) is the partial derivative of kynis with respect to 0^. rjt in Eq. (fT4l) can be 
solved backward given r/i = V(j,-^E, where E = J^^^E^dx, which will be discussed in the 
following. 

Gradient of E with respect to 0i: The computation of V^^-E is not straightforward and 
the Riemannian structure of ODFs has to be incorporated. Let's first compute V^^-E^ at a fixed 
location, x. We consider a variation (pl = (pi + eh of 0i and denote the corresponding variation in 
A as y4% where A = -D^^i and A'' = Dx(t)\ ■ Here, we directly give the expression of deEx\e=o 
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and the reader is referred to Appendix |A] for the full derivation of terms (A) and (B) in the 
following equation. 



deE,U=o (15) 

\^—n > 

\^—r\ I 

A^pteInp04>l^{x) 
Alptcmp04>i^{x) 



n/i I r \ ^ vtGmp"i,v^i ; 
=2^1og^Vtcn.pO0r^W ^targ(x), ^^ |.=0 



2\logAV..o.pO,^-^(.) V'targ(x), "" ' Q^ le=0 



-2\logAVtc„.pO<^-^(x) V'targ(x), ^^ |.=0 



term (A) 



^iogAv.ton,po<^r'(=^) ^''^te'^p ° (<^i) ^(^) 



2\l0gAV>te.np0<^-i(x) V'targ(x), ^ 



AV'tcmpO(/)j^ (x) 



term (B) 



2det0i(x)<^ {{D^(t)i) "^(log^^^^^^^p(^.)Vtarg(0l(x)),V^.(AV'tcmp)>^^^^^p(^),/i 



V 

term (A) 



+ I](div((log^vto.„p(x) V'targl^ila:)), L^>^^^^^_^(^.))e\ /i 



i=l 



term (B) 

T 



where denotes the matrix transpose and e* is a 3 x 1 vector with the zth element as one and 
the rest as zero. (-, ■) A^|)tcrnp{x) is the Fisher-Rao metric defined in Eq. Q. Vs(y4'j/'tcmp) in term 
(A) is the first derivative of the VODF, Ai/jtemp, with respect to x. Since Ai/^tcmp also lies in the 



Riemannian manifold of VODFs, V x{A.il)x.cnip) is a vector with each element being a logarithm 
map of Ai/jtcmp and is defined as 

JK^\ logAv>tcn.p(x) ^V'tcmp(a; + Aei) 
JK^\ logAv>tc„.p(x) ^V'tcmp(a; + Aes) 
\K^\ logAV.tc„.p(x) ^V'tcmp(a; + Aeg) 
where Aei, Ae2 and Aea indicate small variations in three orthonormal directions of M^, respec- 
tively. 



Va; [AV'temp(2;)] 
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12 



In term (B) of Eq. (fT5l) . we define L^. as a 3 x 3 matrix of logarithm maps with its ith column 
written as 

where w* is the ith column of {Dx4>i)^^- Denote 's = (Dj-^i) s. m* is the ith element of vector 



u 



det(D,0i) \D.,<Pr)-^Vs 



V^ 



In sum, V^i-E can be computed by integrating V ^^E^ over the image space and written as 



V<^,i^ 



2 / det(0i(x)) { (Z),.(/.ij 



(16) 



( logAVton.p(a=) V'targ(</'l(a;)) , ^ logAV,temp(xO ^V'tcmp(x + Aci) ) Arp,,^^{x) 
{ logAV,te.np{a;) V'targl^llx)) , ^ log^^^_^(^.) AVtcmp(a; + Aca) )AV>tcn.p(x) 
( ^^^A^tompix) '^targiM^)) , ^ log^^^^^p(^) AV'tcmp(a; + Acg) )AVtcmp(x) 



+ 



> dx . 



div((l0g^^^_^(,.) Vtarg(0l(x)), ^i>^^^_^(,)) 
div((l0g^^,,„,^(,.) V'targ(0l(x)), L2)^^^^^^^^^.^) 
_ div((log^^^^_^^^(,) Vtarg(0l(x)),i:3.)^^^^^^^^^^) J ^ 

We now like to emphasize the difference of this above gradient derivation from our previous 
work [[34ll . The fundamental difference is that in ll34ll . we assume that A does not change under 
the variation (pl and thus, do not consider the variation in A, i.e.. A'' is ignored. Therefore, in 
ll34l . the gradient of E with respect to (pi only incorporates term (A) of Eq. (fT5l) . This term is 
similar to the scalar image matching case and only takes into account image shape difference 
in the volume space. We illustrate this n Figure [2l where we have one template image and two 
target images. Figure |2] (a) shows the template image, where its overall image shape is circular 
and the ODFs at each voxel inside the circle are oriented horizontally. Figure [2] (b) shows the 
first target image, where its overall image shape is an ellipsoid and the ODFs inside its voxels 
are oriented horizontally. Figure |2] (c) shows the second target image, where its overall image 
shape is circular as the template image but the ODFs at each voxel inside the circle are oriented 
at 45°. The results obtained using only term (A) as proposed in [[34ll are shown in Figures |2] (f, g). 
In Figure |2](f), we see that because of the contribution of term (A) in Eq. (fT5l) . the deformation 
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field and its corresponding momentum in the target space point to the direction that enlarges the 
circle to the ellipsoid. However, in Figure [21 (g), we see that term (A) in Eq. (fT5l) is unable to 
account for such deformations as the image shapes are the same, resulting in the deformation 
field being zero. Figures [2] (d, e) show the results using both terms (A) and (B) as proposed in 
this current paper. From Figure |2] (d), we see that the proposed algorithm gives a deformation 
field that enlarges the circle to the ellipsoid, similar to that of Figure |2] (f). More importantly, 
as shown in Figure [2] (e), we see that the deformation that amounts to rotating the ODFs is 
captured by term (B) of Eq. (fTSl) . which is a property that [[34ll does not possess. 

F. Numerical Implementation 

We so far derive J and its gradient VJ{mt) in the continuous setting. In this section, we 
elaborate the numerical implementation of our algorithm under the discrete setting, in particular, 
the numerical computation of W^j^E. 

In discretization of the spatial domain, we first represent the ambient space, Q, using a finite 
number of points on the image grid, f2 = {(xj)^^}. In this setting, we can assume rrit to be the 
sum of Dirac measures such that uit = X]j=i '^i(^) ® ^</>t(xi) such that 



in n 



■^0 .=1 ,=1 



y((/)i(xi),(/)i(a;j))aj(t)], 



where ai{t) is the momentum vector at Xi and time t. In discretization of the spherical domain S^, 
we discretize it into Ns equally distributed gradient directions on the sphere. For each gradient 
direction k, it can be represented as 3D vector with unit length s^ in Cartesian coordinate 
and (rfc, 9k, ^k) in the spherical coordinate. We use a conjugate gradient routine to perform the 
minimization of J with respect to ai{t). We summarize steps required in each iteration during 
the minimization process below: 

1) Use the forward Euler method to compute the trajectory based on the flow equation 

d(l)t{xi) 



dt 



^/cy(0t(xi),0t(xj))aj(t) . (17) 



2) Compute V (i,^(^x.)E in Eq. (fT6l) . which is described in details below. 

3) Solve rjt = ['r]i{t)]fLi in Eq. (fT4l) using the backward Euler integration, where i indices Xj. 

4) Compute the gradient VJ{ai{t)) = 2ai{t) + r]i{t). 
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Fig. 2; The first and second rows respectively illustrate the original HARDI and their enlarged images. Compared to the image 
on panel (a), the image on panel (b) has the same ODFs but a different ellipsoidal image shape, while the image on panel 
(c) shows different ODFs but the same circular image shape. Panels (d) and (e) show the deformations and the corresponding 
momenta, calculated using V^jS in Eq. l ll6t . for mapping the image on panel (a) to panels (b) and (c), respectively. Panels 
(f) and (s) show the deformations and the corresponding momenta, calculated using the gradient in our previous work 1341 . for 
mapping the image on panel (a) to panels (b) and (c), respectively. 
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5) Evaluate J when ai{t) = af'^{t) — eVJ{ai{t)), where e is the adaptive step size determined 
by a golden section search. 
Since steps 1,3 — 5 only involve the spatial information, we follow the numerical computation 
proposed in the previous LDDMM algorithm ||43l . 

We now discuss how to compute V ^^i^Xi)E in Eq. (fT6l ). which involves the VODF interpolation 



in the spherical coordinate for Ai/7temp(a^i) at a fixed Xj and the a/ODF interpolation in the 
image spatial domain for il)x.a.Tg{4>i{^))- To do so, we rewrite Ail)tcmp{xi) as Ai/jtompls/cXj) 
and V'targ(0i(2^i)) as i/'targ(sfe, 0i(a;j)). For the VODF interpolation in the spherical coordinate 
for Ai/?tcmp(a^i) at a fixed Xj, we compute Ai/?tomp(sfc,a:j) according to Eq. (|7]) using angular 



interpolation on §^ based on spherical harmonics. For the yODF interpolation in the image 
spatial domain for 'j/'targ(0i(2^))j we compute Vtarglsfc, (j)i{xi)) under the Riemannian framework 
in §II-AI as 

Vtarg(s.,0l(x,)) =exp^^ (^^^^^^^,^^^ $:^,l0g (^^^^^^^,^^)(Vtarg(s„X,)) , 

where Mi is the neighborhood of Xj, and Wj is the weight of Xj based on the distance between 
01 (xj) and Xj. The exponential maps and logarithm maps can be computed via Eq. (U) and Eq. 
© respectively. Finally, the inner product in Eq. (fT6l) . 

( logAVtc„.p(x) V'targl^lla;)) , ^-— log^^^_^(^.) AVtcmp(x + Ac^) ) AV>tcn.p(x), 

and 

(logAV..o.p(a=)V'targ(0l(x)) , 4>^^,_^(,) 

can be computed using Eq. Q, where Ae^ is the voxel size. 

III. Results 

In this section, we illustrate how LDDMM-ODF performs on both synthetic and children 
brain HARDI data and then compare its performance over registration based on using diffusion 
tensors or fractional anisotropic (FA) scalar images. 
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A. Synthetic Data 

We first illustrate that the HARDI model is useful to align crossing fibers, especially when 
crossing fibers have equal orientation distributions. To do so, we construct two synthetic datasets, 
template and target, where there are two identical fibers perpendicularly crossing each other 
(Figure [3] (a, b)). The orientations of the two crossing fibers differ from the template image 
(Figure [3] (a)) to the target image (Figure [3] (b)). We will compare the performance of LDDMM- 
ODF to the LDDMM algorithm based on DTI (LDDMM-DTI) [gll. We refer the reader to (H 
for detailed mathematical derivation for LDDMM-DTI. 

In the HARDI model, such orientation differences are encoded by the ODFs, while in the 
DTI model, the diffusion tensors of both the template and target data look like disks, where 
the first two eigenvalues being equal and the third eigenvalue being almost zero. Although the 
overall image shapes are the same in both the template and target HARDI data, the LDDMM- 
ODF algorithm is able to characterize the orientation difference of the ODFs between them 
by generating the deformation shown in Figure [3] (d) with the help of term (B) of Eq. (fTSi) . 
LDDMM-DTI fails to find any deformation (Figure [3] (e)) even though the reorientation of the 
tensor is taken into account in the tensor mapping. 

B. HARDI Data of Children Brains 

In this section, we apply our proposed algorithm to real HARDI data. We evaluate the mapping 
accuracy of our LDDMM-ODF algorithm by comparing it with the LDDMM-image mapping 
based on FA (LDDMM-FA) and the LDDMM-DTI mapping based on diffusion tensors using 
the brain datasets of 26 young children (6 years old). All three algorithms are developed under 
the LDDMM framework as given in §II-DI with the exception that the matching functional, E, 
is the least square difference between two image intensities for the image mapping, LDDMM- 
FA, and the Frobenius norm between two tensors for the DTI mapping, LDDMM-DTI. More 
precisely, LDDMM-FA is based on the method developed by MM and LDDMM-DTI is based 
on the method developed by ll44ll . In our implementation however, we optimize the deformation 
with respect to the momentum rather than the velocity (see p2l|). It is important to note that 
all three mapping algorithms used in the following evaluation have the same numerical scheme, 
such that any potential errors due to numerical related issues are avoided and we can make a 
fair comparison. 
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Fig. 3: Comparison between the LDDMM-ODF and LDDMM-DTI algorithms. Panels (a, b) respectively show the template and 
target HARDI and their enlarged images, where the ODF or diffusion tensor at each location contains two crossing fibers with 
equal orientation distribution. Panel (c) illustrates the template HARDI image transformed via the deformation given in panel 
(d), the result of the LDDMM-ODF algorithm. Panel (e) illustrates no deformation found via the LDDMM-DTI algorithm and 
thus the template HARDI image remains. 
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Our image data are acquired using a 3T Siemens Magnetom Trio Tim scanner with a 32- 
channel head coil at the National University of Singapore. Diffusion weighted imaging protocol 
is a single-shot echo-planar sequence with 55 slices of 2.3mm thickness, with no inter-slice 
gaps, imaging matrix 96 x 96, field of view 220 x 220mm^, repetition time=6800ms, echo 
time=89ms, flip angle 90°. 61 diffusion weighted images with b=900s/mm^, 7 baseline (bO) 
images without diffusion weighting are acquired. Notice that the b-value used in our acquisition 
is relatively low when compared to HARDI acquisition where b > lOOOs/mm^ typically. This is 
because the water diffusivity is in general faster in young children's brain than in adults' brain. 
The large b-value could result in significant loss of diffusion signals. In addition, our dataset is 
for the purpose of the comparison between the HARDI and DTI models. Thus, the b-value is 
determined by balancing the needs of both HARDI and DTI acquisition. In the data processing, 
DWIs of each subject are first corrected for motion and eddy current distortions using affine 
transformation to the bO image (where there is no diffusion weighting). We randomly select one 
subject as the template in this study and first align the remaining subjects to this template using 
the affine transformation computed based on the bO images of the subject and the template. 
Then, the DTI is computed using least square fitting [i46| and the FA is calculated from the 
DTI, and the ODF, "^affine transformed > is estimated using the approach proposed in ll25l . We then 
respectively employ the LDDMM-FA, LDDMM-DTI, and LDDMM-ODF algorithms to register 
all subjects to the template. To ensure a fair comparison, we fix the general setting of LDDMM 
with kernel ay = 5 (Eq. ^). For LDDMM-FA and LDDMM-DTI, based on the diffeomorphic 
mappings computed in each case, we apply the diffeomorphic group action defined in Eq. ^ 
to Vaffine transformed to obtain the registered ODFs. 

To evaluate the mapping accuracy for the whole brain, we compute symmetrized KuUback- 
Leibler divergence (sKL) of the ODFs between the deformed subject and the template. The sKL 
has been used as a metric for comparing ODFs in ifTTl and is defined as 

sKL(pi,p2)=/ Pi(s)log^rfs+ / p,{s)log^ds, (18) 

isG§2 P2(s) ise§2 Pi(s) 

for two ODFs pi(-) and P2{-)- Lower sKL indicates that the ODF of the subjects are better 
aligned. Figure |4] illustrates the averaged sKL maps across all 25 subjects when affine, LDDMM- 
FA, LDDMM-DTI, or LDDMM-ODF are appUed. This figure suggests that LDDMM-ODF is 
the best mapping among all studied in this paper as it has the least amount of variation, even 
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though we do not use the sKL metric in LDDMM-ODF. Figure [5] also shows the cumulative 
distributions of sKL across the image space from each mapping. Kolmogorov-Smirnov tests 
on the cumulative distributions also suggest that the LDDMM-ODF significantly reduces sKL 
distance against the other three methods {p < 0.001). 

We now evaluate the mapping accuracy of individual white matter tracts using 1) sKL of the 
ODF between the template's and deformed subject's tract and 2) Dice overlap ratio to quantify the 
percentage of the overlap volumes between the template and deformed subject's tracts. We extract 
three major white matter tracts, including the corpus callosum (CC) and bilateral corticospinal 
tracts (CST-left, CST-right), using probabilistic tractography with the help of Camino [|46l . The 
probabilistic tractography is performed on the q-ball reconstruction using spherical harmonic 
representation up to order 6 with the number of directions for each ODF limited to 3 and the 
maximum allowed turning angle limited to 70°. 

We adopt the anatomical definition of the CC, CST-left and CST-right given in [|47l and define 
three regions of interest (ROI) such that each tract is comprised of all fibers passing through 
these three ROIs. Figure [6] shows the sKL maps for the three tracts, suggesting that, again 
LDDMM-ODF provides the best alignment for the ODFs of these three tracts when compared 
to affine, LDDMM-FA, and LDDMM-DTL Figure |7] shows the average sKL values for the 
CC, CST-left, and CST-right. Moreover, Figure [8] shows the averaged Dice overlap ratios across 
all 25 subjects for the CC and bilateral CST One-sample t-tests shows that LDDMM-ODF 
significantly improves the alignment of local fiber directions for three fiber tracts against the 
other methods in terms of sKL (p < 0.001). In addition, the one-sample t-tests between any 
two mapping algorithms suggest that all the non-linear methods show significant improvement 
against affine in terms of Dice overlap ratio (p < 0.001) for the three tracts, and LDDMM-ODF 
shows significant improvement against LDDMM-FA and LDDMM-DTI {p < 0.001). In the 
comparison between LDDMM-FA and LDDMM-DTI, the only significant difference is found 
in the CST-left (p < 0.05), while no significant differences are found in the CC and CST-right. 

IV. Conclusion 

We present a novel diffeomorphic metric mapping algorithm for aligning HARDI data in 
the setting of large deformations. Our mapping algorithm seeks an optimal diffeomorphic flow 
connecting one HARDI to another in a diffeomorphic metric space and locally reorients ODFs 
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Fig. 4: Panels (a-d) respectively show the maps of mean symmetrized Kullback-Leibler (sKL) divergence of the ODFs between 
the template and the subjects deformed via affine, LDDMM-FA, LDDMM-DTI, and LDDMM-ODF. 
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Fig. 5: sKL Cumulative distributions across the whole brain image and averaged over all 25 subjects are shown in blue for 
affine, cyan for LDDMM-FA, yellow for LDDMM-DTI, and red for LDDMM-ODF, respectively. 
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Fig. 6: Panels (a-h) show the maps of mean symmetrized KuUback-Leibler (sKL) divergence of the ODFs between the template 
and the subjects deformed via affine, LDDMM-FA, LDDMM-DTI, and LDDMM-ODF for the three major white matter tracts 
of the corpus callosum (CC) and bilateral corticospinal tracts (CST-left, CST-right). 



July 26, 2011 



DRAFT 



IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. X, NO. X, XXX 2011 



22 



0.4 

0.35 

0.3 

0.25 

i 0.2 

0.15 

0.1 

0.05 





■affine 

□ LDDMM-FA 

□ LDDMM-DTI 
■ LDDMM-ODHI 




CST-left 



CST-right 



Fig. 7: sKL averaged over all 25 subjects are shown for the corpus callosum (CC) and bilateral corticospinal tracts (CST-left, 
CST-right) when affine (blue), LDDMM-FA (cyan), LDDMM-DTI (yellow), or LDDMM-ODF (red) are applied. 
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Fig. 8: Dive overlap ratios averaged over all 25 subjects deformed by affine (blue), LDDMM-FA (cyan), LDDMM-DTI (yellow), 
or LDDMM-ODF (red). 
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due to the diffeomorphic transformation at each location of the 3D HARDI volume in an 
anatomically consistent manner. We incorporate the Riemannian metric for the similarity of 
ODFs into a variational problem defined under the LDDMM framework. The diffeomorphic 
metric space combined with the Riemannian metric space of ODF provides a natural framework 
for computing the gradient of our mapping functional. We demonstrate the performance of our 
algorithm on synthetic data and real brain HARDI data. This registration approach will facilitate 
atlas generation and group analysis of HARDI for a variety of clinical studies. We are currently 
investigating the effects of our registration algorithm on fiber tractography. 

Appendix A 

Gradient of E^ with respect to 0i 

We now elaborate the derivation of terms (A) and (B) in Eq. (fT5l) . 
Term (A): For the sake of simplicity, we denote term (A) of Eq. (fT5l) as Ea and rewrite 



Given ^Mi±|l^^|,=o = - [{D(j)i)-'h] o 0-i(x), we have 
With a change of variable from x to (pi^{x), we have 

Ea = 2 det (/.i(a;)<^(D^.0i)-^(log^^^_^(^) Vtarg(</'l(a^)), V^.(AVtemp)>^^^_^(^.), h 

Term (B): We denote term (B) of Eq. (fT5l) as Eb and rewrite 






= -2\logAVtc.npo0-i(x)V'targ(a;), ^ 

According to Theorem III. 11 we have 



AV'tcmpO?!>i (x) 



^''0temp(s,x) O0^^(x) 



— ^ ^^ „.^ ,. -1 „ '^ 



\ ||(/^..0i)^'s NIPx</'0 



s 



0r'(x) . 
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Denote s = (D,^,) "'s. Given ^^"'-^^^^"'^"^ 



_ d{Da:<f>i+iD^h)-'^[x) , 



■{DM^^D,h{DM^\ 



we can now compute 



de 



det(D,.0i) Vv^ 



e=0 

''^Vllfil'^) 



v^ 



(D,0l)-^D,.M/^x.0l)-'s 



O01 [X] 



where 



-A'j/)tomp(s,x)trace(^L'^/i(i:'^0i) ^ 
D^h{D^(t)i)~^s,uj - -A'j/'tcmp(s,a;)trace(^L'^./i(L'^.0i)"^j > o0-^i 



X 



det(D,.0i) \z},.0i)-Tv^ 



v^MP 



We now derive the above equation in order to express it in an explicit form of h. Before doing 
so, we first define a 3 x 3 identity matrix as Idsxs = [e\ e^, e^], where e* is a 3 x 1 vector with 
the zth element as one and the rest as zero. Denote (D^0i)^^ = [w^, w^, w"^], where w* is the 
i\h column of (D^</)i)~^. Thus, the trace of Dxh{Dx(j)i)~^ can be written as 

3 

trace (d,/i(D,.0i)-i) =^(D,/iw\e^ 



i=l 



It yields 



M'V'tcmp(s,a;) o0^i(x) 



de 



|e=0 



= J(D,./i(Z},.0i)"'s,n)-^^AV'tcmp(s,x)(z},/iw\e*)io0-i(x) . 
We introduce the following lemma l|48l that leads to a simple expression of Eb- 
Lemma A.l. For smooth vector fields, h, u, w, defined in a bounded open domain in M?, 

diviu^w) , h 



where u" is the ith element of u. 



div(M'^U7) 
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E 



h )o(j)~\x) 



As a consequence, when defining L^ = {D^cjyi) ^sm* — i/l'j/'temp(s,a;)w* , it can be easily 
shown that 

div((log^^^^„^^(,) V'targl^ila;)), Ll) ^^^^^^^^^^) 

With a change of variable from x to (j)^^{x), we finally have 

div((log^^^_^(,) Vtarg(0l(x)), ^i>^^,_^(,)) 
div((log^^^_^(,) Vtarg(0l(x)), ^x>^^,_^(,)) 
div((l0g^^^_^(,.) Vtarg(0l(x)), ^'>^^,_^(,)) J 



EB = 2det(0i(x)) 



M 
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